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Abstract 


A technique is implemented for computing hypersonic aeroheating, shear stress, and other flow 
properties on the windward side of a three-dimensional (3D) blunt body. The technique uses a 
2D/axisymmetric flow solver modified by scale factors for a corresponding equivalent 
axisymmetric body. Examples are given in which a 2D solver is used to calculate the flow at 
selected meridional planes on elliptic paraboloids in reentry flight. The report describes the 
equations and the codes used to convert the body surface parameters into input used to scale the 
2D viscous shock layer equations in the axisymmetric viscous shock layer code. Very good 
agreement is obtained with solutions to finite rate chemistry 3D thin viscous shock layer 
equations for a finite rate catalytic body. 


Introduction 

To assess the design of reentry spacecraft with respect to aeroheating, a number of techniques are 
available, depending on the actual geometry and the needed accuracy. For highest fidelity, 
solutions to the full Navier-Stokes equations would be required, but this is very time-consuming 
both in human labor as well as computer time. It is not practical to use the Navier-Stokes 
solutions for quick assessments of proposed vehicle shapes because of the amount of time it 
takes to define an adequate flow field grid and to run the solution until convergence. The 
purpose of this work is to implement a rapid approximate technique for predicting aeroheating on 
the windward side of three-dimensional (3D) vehicles at angle of attack. The technique involves 
modifying 2D flow equations for an axisymmetric body at zero angle of attack using scaling 
parameters in the equations that accounts for 3D effects in the flow. The method was developed 
originally by Brykina, et al. 1 ' 2 and the factors only depend on body geometry and Reynolds 
number. 

In preliminary stages of design, an approximate geometry and correlations of the heat flux on 
such geometries may be sufficient. Such tools as Detra, Kemp and Riddell 3 , or Fay and Riddell 4 
on the stagnation point of spheres may be adequate, or similarly there are relations for swept 
cylinders. Likewise, flat plate correlations may be used on certain wing shapes or flaps. The 
MINIVER or LANMIN S code can be used for performing approximate heating calculations for 
simple shapes. The INCHES code by Zoby and Simmonds 6 is an axisymmetric technique that 
relies on the Maslen technique for shock shape. A somewhat higher fidelity engineering code 
that also takes into account variable entropy edge conditions is the AEROHEAT code. 7 It uses 
axisymmetric analog technique for 3D boundary layers. There have been some adequate 
attempts to use axisymmetric codes for simple shapes such as the viscous shock layer (VSL) 
method ' and even apply them to complex shapes at angle of attack, by approximating the shape 
as an hyperboloid of revolution. 10 1 1 Higher-order approximations have been developed for 3D 
shapes such as the 3D VSL method and 3D parabolized Navier-Stokes. One method for 
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avoiding using full Navier-Stokes is to solve the 3D Euler (inviscid) flow equations for the shock 
layer flow around a body, then use the results as boundary layer edge conditions for an 
axisymmetric boundary layer solution. 12 ' 13 Even this quasi-3D technique is very labor- and 
computer-intensive. The technique developed in this report is one of intermediate complexity, in 
which 2D axisymmetric methods may be used with modifications to approximate the flow field 
and heat flux on a 3D geometry. This technique yields very good approximate heating that 
includes surface catalytic effects for chemically nonequilibrium flows. Likewise, one also 
obtains other flow properties such as species concentrations near the surface. Although the 
method can be used to scale any set of 2D viscous flow equations, it is applied here to the 
chemically nonequilibrium 2D VSL equations, specifically, the Miner and Lewis code 14 as 
modified to include finite wall catalysis boundary conditions. 15 

This report presents scale factors used to modify the VSL equations written in body-oriented 
coordinates. The method requires computing scaling parameters based on body surface 
derivatives. The code CONVERT computes these scale factors and other geometry parameters 
needed in the VSL3D code*. Paraboloids of elliptical cross section are used here to compare 
with previous work and to illustrate how the method is used. The paraboloid geometry is defined 
by an auxiliary program, PARG. For each meridional plane and angle of attack of interest, it 
computes the coordinates and surface derivatives that the code CONVERT uses. To validate this 
approximate technique we compare its results with 3D results. The Miner and Lewis code solves 
the 2D VSL equations in body-surface and body-normal coordinates (S,r|), whereas the original 
3D body is defined in Cartesian coordinates. Therefore, a third code, INTPX, is used as a post- 
processor. It uses the output coordinates from the VSL3D code to interpolate the original body 
coordinates to obtain corresponding locations in the original coordinate system. This is useful 
for locating the flow properties at locations on the original 3D body. 

Results are obtained for several cases for paraboloids of elliptic cross section and compared with 
3D VSL solutions. 


Analysis 

Body Coordinates and Input Parameters for VSL3D Code 

The modified 2D VSL equations are solved in a streamwise direction along the surface of an 
axisymmetric body that is equivalent to specific meridional planes tp of the 3D body. There is an 
equivalent axisymmetric body (EAB) for each desired meridional plan and angle of attack. The 
meridional plane is defined by a plane through the stagnation point, parallel to the wind velocity. 
The windward symmetry plane represents tp=0 and the leeward symmetry plane corresponds to 

* In this report we denote the axisymmetric viscous shock layer code as VSL and the one modified to handle the 3D 
EAB as VSL3D. 
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<p=180°. The input parameters required by the VSL3D code are the axial distance from the 
stagnation point z n , the distance from the axis to the body r n , the surface distance s, the local 
body curvature in the r n ,z n -plane, K, and the angle 0, which is the local body tangent with respect 
to the velocity vector Vj n f. See Figure 1 for a pictorial definition of these parameters. In 
addition, the modified VSL code requires the scale factor H/H s , which is the ratio of mean 
curvature of the 3D body to that of the EAB. The code CONVERT computes these parameters 
from coordinates and surface derivatives for each meridional plane of the 3D body defined in a 
Cartesian system of coordinates by the equation z=f(x,y). The pitch plane (plane of symmetry) is 
defined by y=0. The geometrical stagnation point is the location on the windward pitch plane at 
which the local surface of the body is normal to the velocity vector. 

For Elliptic Paraboloid 

k=(a/b) 2 Z n 



Figure 1 Coordinate systems for 3D and equivalent axisymmclric bodies. 


The mean curvature of the EAB is given by 

K + r ’cos# 

2 

and the mean curvature of the 3D body at a point on the surface is 
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These expressions are obtained using differential geometry and their specific forms depend on 
the coordinate system used. The scale factor H/H s is applied to the characteristic length or to the 
Reynolds number that appears in the nondimensional form of the flow equations. It is derived in 
Ref. 1 . 


Surface coordinates and surface derivatives in the Cartesian frame are inputs to the code 
CONVERT, along with the corresponding angle of attack and the meridional angle. To calculate 
these parameters, it was necessary to deal with several special cases where certain factors are 
singular or ill conditioned because of the local geometry. It is also necessary to define these 
quantities at the stagnation because it is a singular point. The following analysis gives the 
different forms of the equations used in the program CONVERT, which computes these 
quantities and the equivalent body radius r and streamwise distance s that the VSL code uses. 


At the stagnation point z n = s = r n = 0 and 0 = n/2. At points away from the stagnation point, wc 
have the following equations for the axial distance z n , and the other transformed coordinates x n 
and y n : 


= (z-z 0 ) since + (jc — jc 0 ) cosce 

(3) 

= (z-z 0 ) cosce — (jc — jc 0 ) since 

(4) 

II 

c 

(5) 


where zoand xo are coordinates of the stagnation point in the original Cartesian coordinate 
system. The EAB origin is the stagnation point (x 0 , 0, z 0 ), and its axis z n is directed along the 
free stream velocity. The coordinates of the stagnation point are defined by the equations 


f xo = tana 
z 0 = f(x 0 , 0). 


The equation of the meridional plane passing through the axis z n is 

y = x n tan#] = [(z - z„ ) since + (x- x () )cosa]tan^ . (6) 


The body angle 0 is given by 


sin 6 = 


f, since + cosce 


(7) 


To find the coordinates r n and s we must integrate the body radius and the arc length along the 
surface of the 3D body in a meridional plane. These quantities are given by 



and 


( 8 ) 
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where 


- 

j d? 


— -land 


= cos 6 


Since these integrands are singular at z n = 0, we perform an integration in the direction different 
from z n from z n = 0 to some small distance away, z n] . These integrals r n (z n ) and s(z n ) are then 
performed as 


[dr f dr 
' - = /*-*■ 


for <p * 90°, or 


if tp = 90 c 


if tp # 90°, or 


y fdr , {dr J 

r ’ = h dy+ ^ dz - 


• }f *. 

/.A 


'f rfv Z [ds 

J= A + A (15) 

if tp = 90°. The point at which we make the crossover is a parameter that we have taken 
arbitrarily to be z n i = 0.1 . These expressions are integrated using Simpson’s rule. 

For the region away from the stagnation point where z n > z n i we integrate along the axial 
coordinate, where we use the relations (10) and (1 1). In this zone we integrate with respect to z n 
from z„i to Z„. 

The surface arc length ds and radial distance dr are determined from the following equations 
where we integrate in the x or y direction. When tp is very near tc/ 2 the differential dr is given by 


dr dr „ (f x sin a + cos«)/ v cosa 

—— = — — /, cos or = — = — — = tan Of cos a 

“■>’ “ z n J(/ x cos or- sin or) 2 + /, 2 ) 


and ds is given by 
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ds ds . . 

— = r cosar = f r cos or . 

dy dz/> ^ 


1 + //+A 2 


= /v 


cos a 


(/ t cos a- sin a) ‘ + / v ) cos# 
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When cp * 7t/2 and for points for which z<z n ] 


— = ~p L l(f x +df v )cosor- sin a] = 
ax dz„ 


(/, since + cos «)[(/, +4 T v )cosq - since] 


^(f x cos a- since ) 2 +/ x 2 
= tan6[ (/ t + Jf v ) cos ce - since] 


(18) 


and 


//v //c 

— = —[(/, +df y )co$a- since] = [(/, + #/ v )cosor-sinor] 
dxdz n 

=[(/, +4T y ) cosce - sin cej / cos# 

If <p = k/2, then at the stagnation point [dr n /dy]o = [ds/dy]o = 1. 
If <p * 7t/2, then 

f x . tO + d 0 /„() 


o+/;+/v) 


(/, cosce -since) +/ v 


(19) 


dr n 
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0 

_dx J 
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f 2 

1 xxO 


cos oc + d^f ^, q 
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We can obtain the curvature in terms of the surface derivatives of the 3D body and angle of 
attack a. The body in a Cartesian coordinate frame is defined as z=f(x,y). The body curvature is 
defined from 


K = 


d6 <i(sin6) 


ds dz n 

When (p * k/ 2, we have 

_ (/* +Cv)[/v cosor - sin a(/ v 2 + 1 )] + + df yy )(f x f y sin a + f y cos or) 

(1 + // +/ v 2 ) 3,2 [(/, + #/ v )cosor-sinor] 

where the full derivative along the meridional plane is 

, dy cos «+ /"sin a 

d = — = — tan (p . 

dx 1 - / v sin or tan 

When cp = k/2 equation (22) reduces to 

[f x cos or - sin or(/ v 2 + 1)]+ / vv / v (/, sin or + cos or) 


(21) 


(22) 


(24) 


K = 


/ v cosor(l + /, + / v ) 


2x3/2 


(25) 
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( 26 ) 


The ratio of mean curvatures for the 3D body to an EAB is given by 


H / H s 


where r n is the distance from 


; > + /„(1 + /, ! )- 2 /„/,/, 

(r + r/' cosflXI + f] + /,*)” 

the axis of the EAB. 


The angle of attack is given in the input as a; the meridional angle is (p. Two forms of equation 
for the body curvature k are used, depending on whether the meridional angle <p = jt/2. If 
(p ^ 7t/2, we use the form 

_ (/«0 +^o/^o) 2 COSQr + ^( /^ + ^o/yvo)/nO (/,. 0 sifl « + COS «) 

° “ (1 + cos«(/„„ + 2J 0 /„ 0 + </„„) (2: 

where 


cos a + f Q sin a 

d o = - — - — ; tan (p . 

1 - / 0 sin atari (p 

If cp is about nil, then (22) reduces to the form 

/no cos a + fyyp 2 Ujo sin a + cos a) 
fyyO cosor(l + f? 0 ) m 

The subscript o denotes the stagnation point. 


(28) 


(29) 


The ratio of average curvatures at the stagnation point is 


0 = 


f XX 0 ( 1 + fyn ) + /vv 0 ( 1 + /«0 ) (/«„ COS 2 O' + / ) COS OT 


2v„(l + // 0 ) 


2 \ 3/2 


2k . 


(30) 


After calculating these variables, they are then normalized with respect to the stagnation point 
curvature: 


Zn = Z n Ko (31) 

r n = r n *K 0 * (32) 

S = s*Ko* (33) 

H/H S =(H/H S )V (34) 

K = k* /Ko* (35) 


where the superscript * denotes the dimensional or unsealed variables. The scaled variables, 
along with the body angle 0, are output from the CONVERT code to be used as input in the 
VSL3D code. They represent the axisymmetric body that is equivalent to the actual 3D body 
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along a meridional plane. Appendix A gives a listing of the computer code CONVERT used to 
generate input for the VSL3D code. 

Modifications to the VSL Code 

The VSL equations for an axisymmetric body are given in Ref. 8, but they must be modified 
slightly for correct application of the EAB. From similarity relations 2 it follows that, to obtain 
heat flux, shear stress, and species concentrations on an actual 3D body, we have to solve 2D 
equations for the EAB with variable Reynolds number. Re* = (H s /H)Re. The VSL equations are 
thus modified by the transformation using the 3D body parameters. The only places in the VSL 
code that require modification are in the calculation of Re or the Reynolds number parameter 
e 2 = 1 / Re and the parameter WREF°cR 0 . These parameters, associated with length scaling, are 
multiplied by H/H s at each point along the surface of the EAB. 

Since the VSL3D code has additional inputs due to its new features, a sample input file is 
included in Appendix B. Appendix C contains the UNIX run stream for a series of cases. The 
output from the VSL3D code is found in Fortran unit 23, a sample of which is found in 
Appendix D. The output from the interpolation code is found in Fortran unit 50, a sample of 
which is found in Appendix E. Note that the nomenclature on output from the VSL3D code is 
different from that found in this paper. Table 1 defines the correspondence between the two 
systems of nomenclature. 

Table 1 Symbol Correspondence 


Symbols used in 
this report 

Symbols in VSL3D and 
INTPX code output 

r n 

RS, rs(I) 

z n 

XB, xb(I) 

s 

S, s(I) 

q 

qw 

K 

ck 

H/H s 

HHS 

T 

tauw 

£ 

eps 


Application to an Elliptic Paraboloid 

To test the coding and the accuracy of the methodology, we have implemented a code that 
computes the shape parameters for a paraboloid having an elliptic cross section. This can be 
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expressed in Cartesian coordinates as z — !/2(x“ + ky“), where ->fk =a/b is the ratio of principal 
axes (x-to-y) of the elliptic cross section (or ratio of principal curvatures at the apex of the 
paraboloid). An auxiliary code (PARG) was used to calculate the shape and surface derivatives 
along planes through the wind axis of the paraboloid at angle of attack. The output of PARG is 
then used as the input of CONVERT to generate the geometry data file for VSL3D. PARG 
generates lines of data, given k, the angle of attack a, and the desired meridional plane (p. The 
windward symmetry plane corresponds to cp = 0 and the leeward symmetry plane to <p = 180°. 
The location of the stagnation point is simply the point where the normal to the surface is parallel 
to the wind axis, or where dz/dx = x 0 = tana and y 0 =0. Then z 0 = '/2X 0 2 = >/2tan 2 a. An array of 
coordinates along a meridional plane on the surface of the paraboloid is set up in equal steps in x 
unless cp = 90 . If cp = 90 , then equal steps in y are taken and all x = xo- When the body surface 
angle becomes parallel to the free stream flow (in the code when coscp < 0.01) the line of points 
is terminated because we can only calculate the flow where the surface is presented to the wind. 
When cp ^ o or 1 80° and if a ^ 0 we have the expression 

1 --/l - & tan 2 (Psin 2 or[jc 2 + 2x cot or- (2 + tan 2 or)] 

, • (36) 

k sin a tan (p v ’ 

or if a = 0 then 


y = x tancp. (37) 

If cp=0 or 1 80°, then y=0. 

We then find the array of z values from 

z = f(x,y) = '/2(x 2 + ky 2 ). (38) 

The surface derivatives are determined analytically and we get the expressions f x = x, f y = ky, 

fxx — 1, fxy = 0, and fyy = k. A listing of the code that generates the paraboloid geometry is given 
in Appendix F. 


Results 

To test the code we have applied it to several flight cases for elliptic paraboloids. Table 2 
describes the flight conditions and Table 3 gives the geometrical parameters used in the 
calculations. The surface catalytic recombination coefficients are taken from Ref. 16; and the 
wall temperature is calculated in the code under the assumption of radiation equilibrium. 
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Table 2 Flight Conditions for Test Case for Elliptic Paraboloids 


Altitude, 

Velocity, 

Density, 

Temperature, 

km 

km/s 

kg/m* 

K 

70 

7.25 

5.91x10 s 

200 


Table 3 Geometry Configuration for Test Cases 


Radius at Stag. Point 
in x-z plane, m 

Angle of 
Attack, a 

Cross Section 
Axis Ratio, k 

Meridional 
Plane, tp 

0.5 

0 

0.25 

0 

0.5 

0 

0.25 

45 

0.5 

0 

0.25 

63.4 

0.5 

0 

0.25 

76 

0.5 

0 

0.25 

90 

0.7 

15 

0.4 

0 

0.7 

15 

0.4 

180 

0.7 

30 

0.4 

0 

0.7 

30 

0.4 

180 

0.7 

45 

0.4 

0 

0.7 

45 

0.4 

180 

0.5 

0 

0.4 

0 

0.5 

0 

0.4 

45 

0.5 

0 

0.4 

90 


The first comparison is a flight case for an elliptic paraboloid at angle of attack of zero. This 
elliptic paraboloid has a ratio of principal axes of 0.5, or k=0.25. The method was applied to five 
meridional planes (p where the PARG code was used to compute the geometry and surface 
derivatives for each tp. The code CONVERT transformed these 3D parameters (coordinates and 
surface derivatives) into EAB coordinates parameters z n , r n , s, k, 0, and H/H s . These parameters 
were input into the VSL3D code (the modified axisymmetric VSL code8) and solutions were 
found. Then the streamwise values were interpolated to transform these locations back to the 3D 
body coordinates using a code given in Appendix G. The results for the first comparison are 
shown in Figure 2, where the solid curves are from the axisymmetric analog solution described 
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here and the open circles are solutions of the 3D VSL equations from Ref. 17, solved by an 
implicit finite difference scheme 18 of fourth-order accuracy in normal coordinates and second- 
order accuracy in longitudinal coordinates as implemented by Shcherbak 1 for 3D viscous flows. 
The agreement is quite good, within about 7% or better. 


CM 

E 




Paraboloid k=0.25 R =0.5 m 

O 

h=70 km V =7.25 km/s 

inf 



r 


Figure 2 Heat flux distribution at various meridional angles on an elliptic paraboloid: 
k=0.25, Ro=0.5 m, h=70 km, V«,=7.25 km/s, catalytic rates from Ref. 15, and radiative 
equilibrium wall. Open circles arc 3D viscous shock layer calculations from Ref. 17. 


1 1 




The second comparison is for the same flight case but with a body of larger radius at the 
stagnation point (R o =0.7 m) with elliptic cross section having k=0.4. The angle of attack was 
varied, and the results for the pitch plane (plane of symmetry) are given in Figure 3. These 
results are compared with 3D VSL solutions from Ref. 2. Here the agreement is good, but tends 
to diminish on the windward side far away from the stagnation point. The shear stress for this 
case is given in Figure 4 and is compared with 3D solutions from Ref. 2. 


50 


a 

Present Method 



x’=-[0.5*(x 2 -TAN(a) 2 )*SIN(a)+(x-TAN(a))*COS(a)] 


Figure 3 Heal flux distribution on symmetry plane of elliptic paraboloid at various angles of 
attack with: Ro=0.7 m, k=0.4, h=70 km. V„=7.25 km/s, catalytic rales from Ref. 15, and 
radiative equilibrium wall. Dashed lines are 3D viscous shock layer calculations from Ref. 2. 


As example of surface atom mass fractions, we show results for an elliptic paraboloid having a 
principal axis ratio of 0.707, or k=0.5 at an altitude of 70 km, velocity of 7.25 km/s, and nose 
radius in the pitch plane of symmetry of 0.7 m. Fig. 5 shows the mass fractions of oxygen and 
nitrogen in three meridional planes. Other flow field properties are available for these equivalent 
axisymmetric bodies calculated by the VSL3D code.* 


* The source code for the modified Miner and Lewis viscous shock layer code (VSL3D) is available from the author 
at ES3, NASA Johnson Space Center, Houston, TX 77058. 
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Paraboloid k=0.4 R o =0.7 AoA=15 
h=70 km V jnf =7.25 km/s Scott Catalycity 



Figure 4 Shear stress distribution on elliptic paraboloid at an angle of attack of 15°: R o =0.5 m 
h— 70 km, V_=7.25 km/s, catalytic rates from Ref. 15, and radiative equilibrium wall. 
Dashed lines arc 3D viscous shock layer calculations from Ref. 2. 
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0.25 


Elliptic Paraboloid k=0.4 Ro=0.5 
Vinf=7.25 km/s h=70 km Scott Catalycity 



0 L — — 1 

0 0.5 1 1.5 2 2.5 3 

r„ 

Figure 5 Mass fraction distributions of N and O on meridional planes of an elliptic paraboloid: 
R o =0.5 m. k=0.4, h=70 km. V_=7.25 km/s, catalytic rales from Ref. 15, and radiative equilibrium. 


Conclusions 

An axisymmetric analog technique has been developed and codes to implement it have been 
described that allows one to calculate heat flux distributions on 3D bodies or 2D bodies at an 
angle of attack (thus in a 3D flow field). The technique follows that of Brykina, et al., 1 and 
results from it have been compared with 3D VSL calculations of heat flux on elliptic 
paraboloids. The technique involved modifying the Miner and Lewis 8,15 nonequilibrium, seven- 
species, 2D VSL code. The results indicate that there is very good agreement between the EAB 
results and the 3D calculations. 

The method also can be implemented in any axisymmetric flow code to simulate a 3D body or 
flow by modifying the equations with the scale factor H/H s and using the EAB coordinates. This 
method can be applied to any blunt geometry that can be described by coordinates and surface 
derivatives. Further applications will be the subject of a future publication. 
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Appendix A - Listing of the program CONVERT that computes equivalent axisymmetric 
coordinates and scaling parameter H/H s from input surface geometry data 

c Program Convert 
c 

parameter ( nn= 1 0 0 0 ) 

dimension z (nn) , fx (nn) , fy (nn) , fxx (nn) , fyy (nn) , fxy (nn) 
dimension r (nn) , s (nn ) , ck (nn) , th (nn) , hhs (nn) 

common/body/ xm (nn) , ym (nn) , x (nn) f y(nn) , f s (nn) , f r (nn) , zn (nn) , n 
external fsx, frx, frz, fsz, fsy, fry , fsxm, frxm 
c 

characters na,me 
real*4 k, lzn 

c Input angle of attack alpha, degrees 

c Input angle fi and number of points n for meridional plane 
c Input file for geometry along meridional plane 

c Input small interval lzn near stagnation point for integral on x 
There should be several grid points within the interval lzn. t 
Integrate wrt x near stagnation point and wrt z elsewhere. 
alpha=0 . 
f i = 0 . 

read (31, *)k, alpha, fi 
write (6 , * ) k, alpha, f i 
f id=f i 
j=l 

read (31, *, end=999 , err=888) x ( j ) , y ( j ) , z ( j ) , fx ( j ) , fy ( j ) , 

* fxx ( j ) , fxy ( j ) , fyy ( j ) 
xm( j ) =-x( j ) 
ym(j)=-y(j) 
j=j+l 

go to 1 

write ( 6 , * ) In cgeomn: Error reading unit 31. 7 

continue 
n=j-l 

write (33, 32) ( x ( j ) , y { j ) , z ( j ) , fx ( j ) , fy ( j ) , fxx ( j ) , 

* fxy ( j ) , fyy ( j ) , j=l, 5) 
write (6,*)' no. of points read' ,n 

format ( 8gl6 . 6 ) 
lzn=0 . 1 
c End of input 
c 

alpha=alpha*l .5707963/90 
f i=f i *1 . 5707963/90 

if (abs(fi-l. 5708) .It. 0.001) go to 26 
a=tan ( f i ) 

26 ca=cos (alpha) 
sa=sin (alpha) 
c 

c Stagnation point 

write ( 29 , * ) ' zn ck th' 
zn ( 1 ) =0 . 
s ( 1 ) =0 . 
r ( 1 ) =0 . 

th(l)=l. 5707963 

if (abs(fi-l. 5708) .It. 0.001) then 

ck (1) = ( fxy (1) **2*ca+fyy ( 1 ) **2* ( fx(l) *sa+ca) ) / 

* ( fyy (1) *ca* (1 . +fx ( 1) ** 2 ) ** 1 . 5) 

else 

da= ( ca+fx ( 1 ) *sa) / (1 . -a*fy ( 1 ) *sa) 

c23 4567890123456789 01234567 89 01234567890123 45678901234567890123456789 012 


c 

c 


888 

999 


32 


17 



ck ( 1 ) = ( ( fxx ( 1 ) +a*da*fxy ( 1 ) ) * ( fxx ( 1 ) +da*a*fxy ( 1) ) *ca+a*da* 

* ( fxy ( 1 ) +a*da* fyy ( 1 ) ) *fyy ( 1 ) * ( fx ( 1 ) *sa+ca) ) 

* / ( ( 1 . +fx { 1 ) **2 ) **1 . 5*ca* ( fxx ( 1 ) +2 . *a*da* fxy ( 1 ) +da**2 * 

* a**2*fyy (1) ) ) 
endi f 

hhs (1) =0. 5* (fxx(l) +fyy(l) * ( 1 . +fx ( 1 ) * *2 ) ) / 

* ( ( 1 . +f x ( 1 ) * *2 ) * * 1 . 5*ck{l) ) 
write(29,*)zn(l) , ck ( 1 ) , th(l}*90. /I. 5707963 

c 

c Meridional plane 

c Compute z-distance (zn) , curvature (ck) and angle (th) : 
c 

do 30 j -2 , n 

zn(j)=(z(j)-z(l) ) *ca- (x(j)-x(l))*sa 
i f { abs (fi-1.5708) .It. 0.001) then 

ck ( j ) = ( fxy ( j ) * ( fx ( j ) *ca-fy ( j ) * *2*sa-sa) +fyy ( j ) *fy ( j ) * 

* ( fx ( j ) *sa+ca ) ) / ( f y ( j ) *ca* ( 1 + fx { j ) * *2 + fy ( j ) **2 ) **1.5) 
else 

da= (ca+fx ( j ) *sa) / ( 1 . -a* fy { j ) *sa) 
ck ( j ) = ( ( fxx ( j ) +a*da* fxy ( j ) ) * ( fx ( j ) *ca- fy ( j ) **2*sa-sa) + 

* ( fxy ( j ) +a*da* fyy ( j ) ) * ( fx ( j ) * fy ( j ) *sa+fy ( j ) *ca) ) / 

* ( ( 1 + f x ( j ) * *2 + fy ( j ) **2 ) **1 . 5 * ( (fx(j) +a*da* fy ( j ) ) *ca-sa) ) 
endi f 

th{ j ) =asin ( ( f x ( j ) *sa+ca) /SQRT ( 1+fx (j ) **2+fy(j)**2)) 
write(29,*)zn(j) / ck{j) , th ( j ) *90 . /I . 5707963 
30 continue 

c 

c Compute body radius (r) and surface distance (s) : 

3=1 

do whi le ( zn ( j ) . 1 t . lzn) 

3 = 3+1 
end do 
33=3 
c 

if (abs (fi-1.5708) .It. 0.001) then 
f r ( 1 ) =1 . 
f s ( 1 ) = 1 . 
do 37 j =2 , j j +1 

fr(j ) = ( sa* f x ( j ) +ca) /sqrt ( (fx(j) *ca-sa) **2 + fy(j) **2) *fy(j)*ca 
fs ( j ) =sqrt ( (l+fx(j)**2+fy(j) **2) / ( ( f x ( j ) *ca-sa) **2+fy(j)**2))* 

* fy ( j ) *ca 

c2 3456789012345678901234567890123456789012345678901234567890123456789012 
37 continue 

do 39 j=2, j j 

call simp (fsy,y(j-l) , y ( j ) ,sy) 
call simp ( fry, y(j-l) , y (j ) , ry) 
s ( j ) =sy+s ( j -1 ) 
r ( j ) =ry+r ( j-1) 

39 continue 
else 

if (a . It. -.00000001) then 
d=a/ca 

f r ( j j ) = ( fxx ( 1 ) +d* fxy ( 1) +d* ( fxy ( 1 ) +d* fyy ( 1 ) ) ) / 

* sqrt ( ( fxx ( 1 ) +d* f xy ( 1 ) ) * *2 *ca* *2+ ( fxy ( 1 ) +d* fyy ( 1 ) )**2) 
f s ( j j ) =fr (jj ) 

xm( j j ) =x ( 1 ) 
do 41 j =2 , j j +1 
i=33-3+l 

da= (ca+fx(j)*sa)/ (l.-a*fy(j)*sa) 
fr ( i ) = (sa* fx ( j ) +ca) /sqrt ( { fx ( j ) *ca-sa) **2 + fy(j)**2) 


18 



* * ( (fx( j ) +a*da* fy ( j ) )*ca-sa) 

f s ( i ) =sqrt ( (l+fx(j) * *2 + fy ( j ) * *2 ) / ( ( fx ( j ) *ca-sa) * *2+fy ( j ) * *2 ) ) 

* ( (fx( j ) +a*da*fy ( j ) ) *ca-sa) 
f r ( i ) =- f r ( i ) 

fs (i) — - f s ( i ) 
xm { i ) =x ( j ) 

41 continue 

do 51 i-2 , j j 

call simp { f sxm, xm ( i — 1 ) ,xm(i) ,ss) 
call simp ( frxm, xm ( i — 1 ) , xm(i) , rr) 
s ( i ) =ss+s ( i- 1 ) 
r ( i ) =rr+r ( i — 1 ) 

51 continue 
else 
d=a/ca 

f r ( 1 ) = ( f xx ( 1 ) +d* f xy ( 1 ) +d* ( f xy ( 1 ) +d* fyy ( 1 ) ) ) / 

* sqrt ( (fxx(l) +d*fxy (1) ) * *2 *ca* *2+ ( fxy { 1 ) +d*fyy (1) ) **2) 
f s ( 1 ) = f r ( 1 ) 

do 40 j =2 , j j +1 

da= (ca+fx ( j ) *sa) / (1 . -a*fy ( j ) *sa) 
f r ( j ) = (sa*fx ( j ) +ca) /sqrt { ( fx ( j ) *ca-sa) **2+fy ( j ) * *2 ) 

* * ( ( fx ( j ) +a*da*fy ( j ) ) *ca-sa) 

fs ( j ) =sqrt ( (l + fx( j ) **2 + fy(j ) **2) / ( ( f x ( j ) *ca-sa) * *2 + fy { j ) * *2 ) ) ' 

* ( ( fx ( j ) +a*da*fy ( j ) ) *ca-sa) 

40 continue 

3 4 format ( 4gl6 . 6 ) 
do 50 3=2,23 

call simp (fsx, x( j -1) , x( j ) , ss) 
call simp ( f rx, x ( j -1 ) , x( j) , rr) 

- s ( j ) =ss+s ( j-1) 
r ( j ) =rr+r ( j-1) 

50 continue 
endi f 
endi f 

Do 60 j = j j - 1 , n 

fr ( j ) -tan ( th ( j ) ) 
fs ( j ) = 1 / cos { th( j ) ) 

60 continue 

write (3 2,*) 7 x ( j ) , zn ( j ) , f r ( j ) , f s ( j ) ' 
write (32, 34) (x ( j ) , zn ( j ) # f r ( j ) , f s ( j ) , j =1 , n) 
do 70 j=jj+l,n 

call simp (fsz, zn ( j j ) , zn ( j ) , sss) 
call simp ( frz, zn ( j j ) , zn ( j ) , rrr) 
s ( j ) =s ( j j ) +sss 
r ( j ) =r ( j j ) +rrr 
70 continue 

Compute a ratio of average curvatures of 3D and axisymm. bodies hhs : 
do 80 j=2,n 

h=0 . 5* ( f xx { j ) * (1+fy ( j ) **2) +fyy ( j ) * (l+fx( j ) **2 ) - 

* 2 * fxy ( j ) *fx( j) *fy(j) ) / (l + f x (j) **2 + fy( j) **2) **1.5 

hs-0 . 5* ( ck ( j ) +cos ( th ( j ) )/r(j) ) 
hhs ( j ) =h/hs 
80 continue 

Do scaling for body with stagnation point radius = 1: 
do 9 0 j =2 , n 

zn ( j ) =zn ( j ) *ck(l) 
r ( j ) =r ( j ) *ck (1) 
s ( j ) =s ( j ) *ck(l) 



ck ( j ) -ck ( j ) /ck (1) 
hhs ( j ) =hhs ( j ) *ck ( 1 ) 

90 continue 

zn ( 1 ) =zn ( 1 ) *ck ( 1 ) 
r ( 1 ) =r ( 1 ) *ck ( 1 ) 
s (1) =s (1) *ck(l) 
hhs { 1 ) =hhs { 1 ) *ck ( 1 ) 
ck ( 1 ) = ck ( 1 ) /ck ( 1 ) 
c 

c Create a file for M.-L. code: 
hhs 1=1 . 

c write(40,*) ' jj alpha, fi, n 7 , j j , alpha, f id, n 

write (40,88) ( zn ( j ) , r ( j ) , s ( j ) , ck ( j ) , th ( j ) 

* , hhs ( j ) , j =1 , n) 

c * , hhsl , j =1 , n ) 

88 f ormat ( 6gl 6 . 6 ) 

stop 
end 

function fsx(xx) 
parameter (nn=1000) 

common/body/ xm(nn) , ym(nn) , x (nn) , y (nn) , f s (nn) , f r (nn) , zn (nn) , n 
call intrp3 (xx, x, fs , n, f sxx) 
f sx= f sxx 
return 
end 
c 

function frx(xx) 
parameter (nn=1000) 

common /body/ xm(nn) , ym (nn) , x (nn) , y (nn) , f s (nn) , f r (nn) , zn (nn) , n 
call intrp3 (xx , x, f r , n, f rxx) 
f rx= f rxx 
return 
end 
c 

function fsxm(xx) 
parameter (nn=1000) 

common /body/ xm (nn) , ym (nn) , x (nn) , y (nn) , f s (nn) , f r (nn) , zn (nn) , n 
call intrp3 (xx, x, f s , n, f sxxm) 
f sxm= f sxxm 
return 
end 
c 

function frxm(xx) 
parameter (nn=1000) 

common/ body/ xm (nn) , ym (nn) , x (nn) , y (nn) , f s (nn) , f r (nn) , zn (nn) , n 
call intrp3 (xx , x , f r , n, f rxx) 
f rxm= f rxx 
return 
end 
c 

function fsy(yy) 
parameter (nn=1000) 

common /body/ xm(nn) , ym(nn) , x (nn) , y (nn) , f s (nn) , fr (nn) , zn (nn) , n 
call intrp3 (yy , y , f s , n, f syy ) 
f sy= f syy 
return 
end 
c 

function fry(yy) 
parameter (nn=1000) 
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common/body/ xm(nn) ,ym(nn) , x(nn) ,y(nn) , fs (nn) , fr (nn) , zn (nn) ,n 
call intrp3 (yy , y , fr , n , fryy) 
f ry= f ryy 
return 
end 
c 

function fsz (zz) 
parameter (nn=1000) 

common / body / xm (nn) , ym (nn) , x (nn) , y (nn) , fs (nn) , f r (nn) , zn (nn) , n 
call intrp3 (zz, zn, fs,n, fszz) 
f sz= f szz 
return 
end 
c 

function frz(zz) 
parameter (nn=1000 ) 

common /body/ xm(nn) ,ym(nn) ,x(nn) ,y(nn) , fs(nn) , fr(nn) ,zn(nn) ,n 
call intrp3 (zz, zn, fr,n, frzz) 
f rz~ f rzz 
return 
end 
c 

SUBROUTINE INTRP3 (XX, X, Y,NPNTS, YY) 

C 

C SUBROUTINE INTRP3 SETS UP THE CALLING ARGUMENT FOR 

C SUBROUTINE INTER3 

C 

C SUBROUTINE INTRP3 CALLS SUBROUTINE INTER3 

C 

C SUBROUTINE INTRP3 IS CALLED BY MAIN. 

C 

C YY IS THE VALUE RETURNED FROM ARRAY Y 

C WHICH CORRESPONDS TO THE VALUE XX IN ARRAY X 

C 

DIMENSION X(NPNTS), Y(NPNTS) 

C 

DATA SMALLT / 1.0D-6 / 

FAC=1 . OdO+SMALLT 
JC=0 

10 JC=JC+1 

IF (XX.GT.X(JC) *FAC) GO TO 10 

IF (JC.LT.2) JC=2 

IF (JC.GT. (NPNTS-1) ) JC=NPNTS-1 

CALL INTER3 ( XX , X ( JC-1 ) , X ( JC ) , X ( JC + 1 ) , Y ( JC-1 ) , Y ( JC) , Y ( JC+1 ) , YY ) 

RETURN 

END 

SUBROUTINE INTER3 ( X , XI , X2 , X3 , FI , F2 , F3 F) 

C 

C SUBROUTINE INTER3 INTERPOLATES FOR THE VALUE F CORRESPONDING TO 

C POINT X USING 3 POINT LAG RANG I AN INTERPOLATION 

C 

C SUBROUTINE INTER3 IS CALLED BY SUBROUTINES INTRP3 , HCP AND HCPA 

c ' 

c ASSUMES XI .LE. X . LE . X3 . 

C 

C WRITE ( 0 , * ) ' INTER3 : ENTRY ' 

C 

AN1 = (X-X2) * (X-X3) 

AN2= (X-Xl) * (X-X3 ) 

AN3 = (X-Xl) * ( X-X2 ) 

DN1= (X1-X2) * ( X1-X3 ) 
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DN2= (X2-X1) * (X2-X3) 

DN3 = { X3 -XI ) * (X3-X2) 

CN1=AN1 /DN1 
CN2=AN2 /DN2 
CN3 = AN3 /DN3 

F=CN1 *F1+CN2 *F2+CN3 *F3 
C WRITE (0,*) ' INTER3 : RETURN' 

RETURN 
END 

c integral by simpson's rule 
c 

subroutine simp ( f unc , a , b, s ) 
parameter ( eps=l . e-2 , jmax=2 0 ) 
os t = - 1 . e3 0 
os = - 1 . e3 0 
do 11 j = 1 , jmax 

call trapzd ( func , a , b, st , j ) 
s= ( 4 . * st-ost ) / 3 . 
i f (abs ( s-os) . 1 t . eps*abs (os ) ) 
os = s 
ost =st 
11 continue 

pause 'too many steps' 
end 
c 

subroutine trapzd ( func , a , b, s # n 
external func 
i f (n . eq . 1 } then 

s = 0 . 5 * (b-a) * ( func ( a ) + f unc ( b) 
i t= 1 
else 

tnm- i t 

del= ( b-a ) /tnm 
x=a + 0 . 5 *del 
sum-0 . 

do 11 j =1 , it 

sum=sum+f unc (x) 
x=x+del 

11 continue 

s=0.5* ( s+ (b-a) * sum/ tnm) 

i t=2 * i t 

endi f 

return 

end 


return 

) 

) 
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Appendix B - Input data file for the modified Miner and Lewis axisymmetric viscous shock layer 
code, VIS3D 


Z-7 0 . km 
$ input 
AL«T= 7 0.00, 
BRAD-1 . 640, 
DS= . 03 , 

I END-2 00 , 
SEND- 3 .0, 
RINF-1 . 14e-7 
TB-2600 . # 
TINF-3 6 0 . , 
UINF-23 7 8 6 . , 
XLE-1 . , 
IGEOM-3 , 


Ro-0.5 m Brykina aoa-0 VINF-7.25 KM/SEC Finite CAT Scott 


! For info and labeling only 

• Nose radius in feet (in pitch plane) 

! Initial streamwise step size 
■ Maximum number of streamwise steps 
! Final S value in nose radius 

• Free steam density, slugs/cuft 
! Wall temperature, R 
! Free stream temperature, R 
! Free stream velocity, ft/s 
! Lewis number 

! Geometry flag. 3-geometry read in from unit 4 
NAN-1 , NDATA-1 , NS- 6 , NSI-6 , I VSL standard inputs, eg. No. of species 
CAT-4, ! Catalycity flag. 4 = compiled in value for Scott's RCG value 

DSN-1 ., RAT-1 . 1 , DSNN- . 004 , ! Streamwise clustering parameters 

THINI-1., ! Shock layer option. 1-Thin Shock layer equations 0-Full VSL 

PRNTCI-0 . , NTSH-2 0 , XKETA-1 . 08 , KEND-1 , KTWAL-0 , ! VSL standard input options 

Flag for Tw. 0-Const, input value. 1-Radiation equilibrium 


ITSW-1 , 
$end 


. 1 


.233 
. 133 


0 . 


.1 


0 . 


.767 
. 677 


Free stream mass fractions 
Initial wall mass f racs . 


Appendix C - Sample output from VSL3D code Fortran Unit 23 


k 

i 

S 

XB 

RS 

ck 

HHS 

eps 

qw 

tauw 

1 

1 

0 

0 

0 

1 

1 

1.98E-02 

50.5 

0 

1 

2 

0.033 

0.0005 

0.033 

1 

1 

1.98E-02 

49.38 

6.572 

1 

3 

0.0693 

0.0024 

0.0692 

1 

1 

1 .98 E-02 

49.38 

13.97 

1 

4 

0.1092 

0.006 

0.109 

1 

1 

1.98E-02 

50.13 

22.82 


cf 

C(O) 

C(02) 

C(NO) 

C(N) 

C(NO+) 

C(N2) 

0.00E+00 

0.1798 

4.46E-02 

1.89E-02 

1 .37E-04 

9.94E-1 1 

0.7565 

8.76E-04 

0.1799 

4.46E-02 

1.88E-02 

1.37E-04 

9.92E-1 1 

0.7565 

1.86E-03 

0.1784 

4.61 E-02 

1.88E-02 

1.29E-04 

9.92E-1 1 

0.7565 

3.04E-03 

0.1776 

4.68E-02 

1.90E-02 

1.25E-04 

9.93E-1 1 

0.7565 


Appendix D - Sample output from INTPX code Fortran Unit 50 


ns=26 

nb=201 

ni=201 




s(i) 

xb(i) 

rs(i) 

X 

y 

z 

0 

0 

0 

0 

0 

0 

0.033 

0.0005 

0.033 

0 

0.132 

0.0022 

0.0693 

0.0024 

0.0692 

0 

0.277 

0.0096 

0.1092 

0.0059 

0.109 

0 

0.4359 

0.0238 

0.1532 

0.0116 

0.1526 

0 

0.6104 

0.0466 

0.2015 

0.02 

0.2001 

0 

0.8007 

0.0801 

0.2546 

0.0317 

0.252 

0 

1 .0078 

0.127 

0.3131 

0.0475 

0.3083 

0 

1.2331 

0.1901 

0.3774 

0.0681 

0.3692 

0 

1.4766 

0.2726 

0.4481 

0.0945 

0.4348 

0 

1 .7386 

0.3778 


r 

zn 

qw 

tau 

cf 

cO 

0 

0 

48.7 

0 

0 

0.192 

0.132 

0.0005 

48.42 

3.134 

2.03E-03 

0.1921 

0.277 

0.0024 

48.46 

6.686 

4.33E-03 

0.192 

0.4359 0.0059 

48.51 

10.53 

6.82E-03 

0.1919 

0.6104 0.0116 

47.87 

14.55 

9.42E-03 

0.192 

0.8007 0.02 

46.81 

18.67 

1.21 E-02 

0.1926 

1.0078 0.0317 

45.32 

22.79 

1.48E-02 

0.1937 

1.2331 

0.0475 

43.42 

26.73 

1.73E-02 

0.1952 

1.4766 0.0681 

41.1 

30.34 

1 97E-02 

0.1972 

1.7386 0.0945 

38.41 

33.43 

2.17E-02 

0.1998 


c02 cNO cN cNO+ cN2 

3.69E-02 1.06E-02 6.62E-02 9.83E-11 0 6943 

3.69E-02 1.06E-02 6.62E-02 9.81 E-1 1 0.6943 
3.69E-02 1.06E-02 6 61 E-02 9.80E-11 0 6944 
3.71 E-02 1.06E-02 6.62E-02 9 82E-11 0.6943 
3 69E-02 1.06E-02 6.64E-02 9.81 E-1 1 0.6941 
3.64E-02 1.05E-02 6.67E-02 9.80E-11 0.6939 

3.55E-02 1.02E-02 6.73E-02 9.79E-11 0 6933 

3.41 E-02 9.78E-03 6.84E-02 9.77E-11 0.6925 

3.42E-02 9.23E-03 7.01 E-02 9.74E-11 0.6911 
3.02E-02 8.53E-03 7.24E-02 9.71 E-1 1 0 6891 
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Appendix E - UNIX run stream for computing a set of cases for a paraboloid. 

parg 

cp fort. 31 f31. 25-0-0 
cgeomn 

cp fort. 40 fort. 4 
cp fort. 40 f40. 25-0-0 
vsl3d<in-Sl . 6tre>ou .25-0-0sl.6tre 
cp fort. 23 f23.25-0-0sl.6tre 
intpx 

cp fort. 50 f 50 . 25-0-0sl . 6tre 
parg 

cp fort. 31 f31. 25-0-45 
cgeomn 

cp fort. 40 f ort . 4 
cp fort. 40 f40. 25-0-45 
vsl3d<in-Sl . 6tre>ou .25-0-45sl . 6tre 
cp fort. 23 f23.25-0-45sl.6tre 
intpx 

cp fort . 50 f 5 0 . 2 5 - 0 -45s 1 . 6 tre 
parg 

cp fort. 31 f 3 1 . 2 5 - 0 - 63 . 4 
cgeomn 

cp fort. 40 f ort . 4 
cp fort. 40 f 4 0 . 25 - 0 - 63 . 4 
vsl3d<in-Sl . 6 tre>ou .25-0-63 . 4sl . 6tre 
cp fort. 23 f23 .25-0-63 . 4sl . 6tre 
intpx 

cp fort. 50 f50. 25-0-63 .4sl. 6tre 
parg 

cp fort. 31 f31. 25-0-76 
cgeomn 

cp fort. 40 fort. 4 
cp fort. 40 f40. 25-0-76 
vsl3d<in-Sl . 6tre>ou.25-0-76sl . 6tre 
cp fort. 23 f 2 3 . 2 5- 0 - 7 6 si . 6 tre 
intpx 

cp fort. 50 f 50 . 2 5- 0-76sl . 6 tre 
parg 

cp fort. 31 f31. 25-0-90 
cgeomn 

cp fort. 40 fort. 4 

cp fort. 40 f40. 25-0-90 

vs 13 d< in- SI . 6tre>ou.25-0-90sl . 6 tre 

cp fort. 23 f 23 . 2 5-0-90sl . 6tre 

intpx 

cp fort. 50 f50.25-0-90sl.6tre 
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Appendix F - Listing of the program PARG that generates geometric data for an elliptic paraboloid 

C Program 3D paraboloid geometry along meridional plane 
c 

parameter (nn=1000 ) 

dimension y (nn) , z (nn) , fx (nn) , fy (nn) , fxx (nn) , fyy (nn) , fxy (nn) 
dimension x(nn) 
c 

real*4 k 

c Input ratio of main curvatures at a stagnation point k 
c Input angle of attack alpha, degrees 

c Input angle fi, number of points n and step dx for meridional plane 
c 

k=0 . 4 
alpha=3 0 . 
fi = 0. 
n=2 01 

write (6,*)' Input k, alpha, fi' 
read (5 f * ) k, alpha, f i 
write ( 6 , * ) ' Input file number 7 
c read ( 5 , *) ifile 

c wri te ( i file ,*) k , alpha , fi 

write (31, *)k, alpha, fi 
dx= 0 . 04 

if ( f i . gt .90.) dx- -dx 
c End of input 
c 

alpha=alpha* 1 .5707963/90 
f i=f i*l . 5707963/90 

i f ( abs (fi-1.5708) .It. 0.001) go to 17 
a= tan { f i ) 

17 ca=cos (alpha) 
sa-sin (alpha) 
c 

wr i te { 6 , * ) 7 ca sa',ca,sa 
c Stagnation point 
x ( 1 ) =sa/ca 
V (1) =0 

z(l)=0.5*x(l)**2 
fx(l)=x(l) 
fy (1) =0 

fxy ( 1) =0 
fxx(l) =1 
fyy (1) =k 
c 

c Meridional plane 
do 21 j=2, n 

cosf i= ( sa*fx ( j -1 ) +ca) /SQRT (1+fx (j-1) **2+fy ( j-1) **2) 
if (cosf i . It . 0 . 01 ) go to 22 

if (abs (fi-1.5708) . It . 0 . 001 ) then 
y ( j ) =y ( j-i) +dx 
X ( j ) =x ( 1 ) 
else 

x ( j ) =x ( j -1 ) +dx 
fx ( j ) =x ( j ) 

if (abs (a) . It . 0 . 000001 ) then 
y(j)=o. 
else 

i f ( alpha . ne . 0 . ) then 

sss= (l-k*a*a*sa* (sa*x( j ) **2 
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* + 2 *ca*x ( j ) -2 * sa-sa* *3 /ca* *2 ) ) 

c write (34,*)j # x(j) , sss 

if(sss.le.O) then 

c write ( 6 , * ) ' In parg : sss is less than 0.' 

go to 777 
endi f 

y ( j ) = ( 1-SQRT (l-k*a*a*sa* (sa*x( j ) **2 + 

* 2 *ca*x ( j ) -2*sa-sa**3/ca**2) ))/(a*k*sa) 
else 

y ( j ) =a*x ( j ) 
endi f 
endif 
endi f 

z ( j ) =0 . 5* (x ( j ) * *2 +k*y { j ) **2) 

f xx ( j ) — 1 

fy ( j ) =k*y { j ) 

fyy ( j ) =k 

fxy ( j ) =0 

j j = j 

777 continue 

21 continue 

22 nfi=jj 

c 

c write (if ile, 32) (x(j),y(j),z(j),fx(j),fy(j), fxx ( j ) , 

write (31 # 32) (x(j),y(j),z(j),fx(j),fy(j), fxx ( j ) , 

* fxy ( j ) , fyy ( j ) . j=l,nfi) 

32 f or mat { 8gl 6 . 6 ) 

stop 

end 


26 



Appendix G - Listing of program I N I PX that interpolates geometrical data by streamwise distances 
calculated in the viscous shock layer code 

c Program INTPX 

dimension x ( 3 0 0 ) , y ( 300 ) , z ( 3 00 ) , rn ( 300 ) , rs ( 3 00 ) , xb ( 3 00 ) , s ( 3 0 0 ) 

dimension zn(300),sn(300),qw(300),ta(300),c(300,5),cf(300) 

character * 8 a 

c Arrays x , y , z , zn , rn , sn have indices that correspond to the same input 
c body points. 

c Arrays s, xb, rs are output from VSL at the same body point, but not same as 
c for the input above. 

c This code interpolates the input arrays to find values corresponding to the 
c output points. 

1 format ( 15x, f 10 . 0 , 2f 15 . 0 ) 
ni = 0 

do 2 i = l , 300 

read (4,3, end=999 , err=888 ) zn ( i ) ,rn(i) ,sn(i) ,ck,th,hss 
ni=ni+l 

2 continue 

888 write ( 6 , * ) In intpx Error detected in read of unit 4' 

999 continue 

read (31, *)k, alpha, fi 

3 format ( 8gl 6.6) 
nb=0 

do 22 1=1,300 

read (31, * , end=9 9 7,err=887)x(i) ,y(i) ,z(i) , f x , fy , f xx , f xy , f yy 
nb=nb+l 
22 continue 

887 write ( 6 , * ) In intpx Error detected in read of unit 31 ' 

997 continue 
ns = 0 

read ( 23 , 10 ) a 
1 0 format ( a8 ) 

do 4 i = l , 3 00 

read ( 23 , 550 , end=998 , err=886)k,ii,s(i) , xb ( i ) , rs ( i ) , ck,hhs,eps 

* ,qw(i) ,ta(i) ,cf (i) , (c(i, j) , j=l,6) 

ns=ns+l 

4 continue 

886 write ( 6 , * ) In intpx Error detected in read of unit 23' 

998 continue 

550 format (2i4,3fl0.4,13gl2.4) 

write ( 50 , * ) 7 ns=',ns,' nb=',nb,' ni= / / ni 

write (50, * ) ' s ( i ) xb(i) rs(i) x y z r zn qw tau cf 

* cO c02 cNO cN cNO+ cN2 
if (ni . gt .nb) ni=nb 

do 20 i=l,ns 

call intrp3 (s (i) ,sn,x,ni,xx) 
call intrp3 ( s ( i ) , sn, y , ni , yy) 
call intrp3 (s (i) ,sn,z,ni,zz) 
call intrp3 (s (i ) ,sn,r,ni,rr) 
call intrp3 (s (i) , sn, zn,ni, zzn) 
rr=sqrt (xx* *2+yy* *2 ) 

write ( 50 , 6 0 ) s ( i ) , xb ( i ) , rs ( i ) , xx , yy , zz,rr, zzn, qw ( i ) , ta ( i ) , 

* cf(i), (c (i, j ) , j = l, 6) 

20 continue 

60 format (3fl0.4,5fl0.4,9gl2.4) 
stop 
end 

SUBROUTINE INTRP3 ( XX , X , Y , NPNTS , YY ) 
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c 

C SUBROUTINE INTRP3 SETS UP THE CALLING ARGUMENT FOR 

C SUBROUTINE INTER3 

C 

C SUBROUTINE INTRP3 CALLS SUBROUTINE INTER3 . 

C 

C SUBROUTINE INTRP3 IS CALLED BY MAIN. 

C 

C YY IS THE VALUE RETURNED FROM ARRAY Y 

C WHICH CORRESPONDS TO THE VALUE XX IN ARRAY X 

C 

c 

DIMENSION X(NPNTS) , Y(NPNTS) 

C 

DATA SMALLT / 1.0D-6 / 

FAC=1 . 0 dO + SMALLT 
JC=0 

10 JC=JC+1 

IF (XX .GT.X ( JC) *FAC) GO TO 10 

IF (JC.LT.2) JC=2 

IF (JC.GT. (NPNTS-1 ) ) JC=NPNTS-1 

CALL INTER3 (XX, X ( JC- 1 ) , X ( JC) , X ( JC+1 ) , Y ( JC-1 ) . Y ( JC) , Y ( JC+1 ) , YY) 

RETURN 

END 

SUBROUTINE INTER3 ( X , XI , X2 , X3 , FI , F2 , F3 , F) 

C 

C SUBROUTINE INTER3 INTERPOLATES FOR THE VALUE F CORRESPONDING TO 

C POINT X USING 3 POINT LAGRANGIAN INTERPOLATION. 

C 

C SUBROUTINE INTER3 IS CALLED BY SUBROUTINES INTRP3 , HCP, AND HCPA. 

C 

C ASSUMES XI .LE. X . LE . X3 . 

C 

C WRITE (0,*) ' INTER3 : ENTRY' 

C 

AN1 = ( X-X2 ) * (X-X3 ) 

AN2= (X-Xl ) * (X-X3 ) 

AN3= (X~X1 ) * (X-X2 ) 

DN1= (X1-X2 ) * (X1-X3) 

DN2= (X2-X1) * (X2-X3 ) 

DN3 = (X3-X1) * (X3-X2 ) 

CN1 =AN1 /DN1 
CN2=AN2/DN2 
CN3 = AN3 / DN3 

F=CN1 *F1+CN2 *F2+CN3 *F3 
C WRITE (0,*) ' INTER3 : RETURN' 

RETURN 

END 
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